# Thursday, 25.May.2006
# Generation of a world map using colors generated by distance function

# Open libraries
library(maptools)
library(foreign)
library(RColorBrewer)


#### for Figure 4 ####
####
dd<-c("C:/XunCao/Research/W_H_backup/MCMCoutputs/MID/dpt/new_latent/Post_Analysis/2000fromsdsc")
setwd(dd)
burn<-c(1000,
        1000,
        2000,
        3000,
        2000,
        2000,
        2000,
        1000,
        2000,
        4000,
        3000)

i<-1 # this is for 2000. change to 9 we have 1990.
Year<-c(seq(1950,2000, by=5)) 
dat<-dget(paste("hyb_", Year[i], "/data", Year[i], sep=""))
Yd<-dat[[1]]
cnames<-rownames(Yd) # gets country names from Yd matrix
n<-length(cnames) # sets n

Za<-read.table(paste("hyb_", Year[i], "/U", sep=""))
dim(Za)
Za<-Za[-c(1:(burn[i]*n)),]

Zb<-read.table(paste("hyb_", Year[i], "/V", sep=""))
dim(Zb)
Zb<-Zb[-c(1:(burn[i]*n)),]



# come up with the color scheme:
dst<-dat[[2]][,,4]
diag(dst)<-0
xy<-cmdscale(dst,k=3)
r<-rank(xy[,1])/dim(dst)[1]
g<-rank(xy[,2])/dim(dst)[1]
b<-rank(xy[,3])/dim(dst)[1]
cc.all<-(cbind(rownames(dst),rgb(r,g,b)))




ZTZ<-matrix(0,n,n)
for(ns in 1:(dim(Za)[1]/n)){
Zns<-as.matrix(Za[ (1+n*(ns-1) ) : (n*ns) ,]) 
ZTZ<-ZTZ+ (Zns)%*%t(Zns)
                 }
ZTZ<-ZTZ/(dim(Za)[1]/n)

ZTZ.v<-matrix(0,n,n)
for(ns in 1:(dim(Zb)[1]/n)){
Zns.v<-as.matrix(Zb[ (1+n*(ns-1) ) : (n*ns) ,]) 
ZTZ.v<-ZTZ.v+ (Zns.v)%*%t(Zns.v)
                 }
ZTZ.v<-ZTZ.v/(dim(Zb)[1]/n)

Z.m0<-sqrt(diag(ZTZ))
ZtZ.n<-ZTZ/( Z.m0%*%t(Z.m0) ) ;  diag(ZtZ.n)<-1
theta<-acos(ZtZ.n)
diag(theta)<-0
xy0<-cmdscale(theta,k=2)  # scales them to two dimensions 
mdwpos<-cmdscale(theta,k=1)
rownames(mdwpos)<-cnames

Z.m0.v<-sqrt(diag(ZTZ.v))
ZtZ.n.v<-ZTZ.v/( Z.m0.v%*%t(Z.m0.v) ) ;  diag(ZtZ.n.v)<-1
theta.v<-acos(ZtZ.n.v)
diag(theta.v)<-0
xy0.v<-cmdscale(theta.v,k=2)  # scales them to two dimensions 
mdwpos.v<-cmdscale(theta.v,k=1)
rownames(mdwpos.v)<-cnames

#postscript("hoffwardjprfigure2.ps",height=3,width=6,horizontal=F,family="Times")
par(mfrow=c(1,2))
par(mar=c(1.5,1.5,1,1))
par(cex=.8)
par(mgp=c(1.75,.75,0))
# plot(xy0*1.4,type="n",xlab="",ylab="",pch=19,col=cc.all[,2],xlim=c(-1.7,1.7),ylim=c(-1,1),main="",las=1,cex=Z.m0*.9) 
plot(xy0*1.4, xlab="",ylab="",pch=19,col=cc.all[,2],las=1,cex=Z.m0*1.2, main="", xlim=c(-2.3,2.3),ylim=c(-1.3,1.3)) 
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
# text(xy0[Z.m0>.5,],cnames[Z.m0>.5], cex=Z.m0*1.5, col=cc.all[,2][Z.m0>.5]) #Z.m0[Z.m0>3.5]*.3
text(xy0[Z.m0>.5,]*1.4,cnames[Z.m0>.5])#, cex=Z.m0*1.5, col=cc.all[,2][Z.m0>.5]) #Z.m0[Z.m0>3.5]*.3

#plot(xy0.v[,1]*2.4,xy0.v[,2]*2.4,type="n",xlab="",ylab="",pch=19,col=cc.all[,2],xlim=c(-1.7,1.7),ylim=c(-1,1),main="",las=1,cex=Z.m0.v*1.3) 
plot(xy0.v[,1]*3.8,xy0.v[,2]*3.8, xlab="",ylab="",pch=19,col=cc.all[,2],xlim=c(-3,3),ylim=c(-2,2),main="",las=1,cex=Z.m0.v*2.3) 
abline(h=0,lty=2,col="gray")
abline(v=0,lty=2,col="gray")
#text(xy0.v[Z.m0.v>.5,]*2.4,cnames[Z.m0.v>.5], cex=Z.m0.v*1.5, col=cc.all[,2][Z.m0.v>.5]) #Z.m0[Z.m0>3.5]*.3
text(xy0.v[Z.m0.v>.5,]*3.8,cnames[Z.m0.v>.5])#, cex=Z.m0.v*1.5, col=cc.all[,2][Z.m0.v>.5]) #Z.m0[Z.m0>3.5]*.3
#text(xy0.v,cnames, cex=Z.m0.v*.3, col=cc.all[,2]) 
# dev.off()


### Then the colored map: 
dp<-c("C:/XunCao/Research/W_H_backup/MCMCoutputs/MID/dpt/Poisson/Pois_Madrid")
setwd(dp)

d0<-read.csv("example1.csv",header=T,sep=",",as.is=TRUE) #base data for map
d1<-read.csv("countrywithcolors.csv",header=T,sep=",",as.is=TRUE,colClasses=c("character")) #geocolors
# d10<-d10[order(d10$idnumber),]      # circumvent merge's reordering
# rownames(d10)<-d10$idnumber       # circumvent merge's reordering 

# read shape file from ESRI
w.shp <- read.shape("robinson.shp") # use equal area projection (Robinson)
w.shp$att.data

Tla<-rep(NA, dim(w.shp$att.data)[1])
for (i in 1:dim(w.shp$att.data)[1]){if (is.element(w.shp$att.data$CNTRY_NAME[i], d0$country.name)){Tla[i]<-d0$tla[d0$country.name==w.shp$att.data$CNTRY_NAME[i]]}}
#color<-rep(rgb(1,1,1), dim(d10)[1])
#for (i in 1:length(Tla)){if (is.element(Tla[i], cc.all[,1])){color[i]<-cc.all[cc.all[,1]==Tla[i],2]}}
#cbind(color,Tla)[1:100,]
#cc.all[1:100, ]

# transform polygons to a map   
wmap <- Map2poly1(w.shp,region.id = as.character(w.shp$att.data$FID))
wmap$plotOrder

par(mfrow=c(1,1))
plot(wmap, #col=color,# forcefill =FALSE,
     xaxt="n",yaxt="n",lwd=.000000000125,
     las=1,ylab="",main="",xlab="") 
     #plot world map with only boundaries

# postscript("hoffwardjprfigure1d.ps",height=5,width=5,horizontal=F,family="Times")
plot(wmap,col=d1$colors,xaxt="n",yaxt="n",lwd=0.0000000001) 
# dev.off() # close graphics device





#### For Figure 6 ####
### let's use the UDV from 10 year intervals:

UTV<-dget("C:/XunCao/Research/W_H_backup/MCMCoutputs/MID/dpt/Poisson/Pois_Madrid/Pois_binaryk0_4_Y/Resultsk3/M.pm.K3.p0.5")
sM<-svd(UTV)

Z<-dget("C:/XunCao/Research/W_H_backup/MCMCoutputs/MID/dpt/Poisson/Pois_Madrid/Pois_binaryk0_4_Y/data9000")[[1]] ; diag(Z)<-NA
m<-dim(Z)[1] ; n<-dim(Z)[2]
K<-2

# postscript("C:\\Documents and Settings\\mward\\Desktop\\rsmdw dp project\\Slides\\doublering2000.ps",horizontal=F,height=6,width=6)
par(mfrow=c(1,1),mar=c(1,1,1,1),oma=c(1,1,1,1) )

act<- (1:n)[  apply(Z,2,sum,na.rm=T)>0  | apply(Z,1,sum,na.rm=T)>0 ] 
Zact<-Z[act,act] 
nact<-dim(Zact)[1]
actors<-rownames(Zact)
u<-sM$u[act,1:2]
v<-sM$v[act,1:2]

ssum<-apply(Zact,1,sum,na.rm=T)
rsum<-apply(Zact,2,sum,na.rm=T)
rs<-sqrt( apply(u^2,1,sum) )
rr<-sqrt( apply(v^2,1,sum) )

xys<-diag(1/rs)%*%u*.6
xyr<-diag(1/rr)%*%v

##

##

plot(xyr,type="n",xaxt="n",yaxt="n",xlab="",ylab="")#,bty="n")
for(i in  (1:nact)[ssum>0] ) {
for(j in ((1:nact)[-i]) [Zact[i,-i]>0 ]  ) 
{ segments(xys[i,1],xys[i,2], xyr[j,1],xyr[j,2],#lwd=Zact[i,j],
           col="lightblue")
                                      } }

text(xys[ssum>0,1],xys[ssum>0,2],actors[ssum>0],cex=2.7*(rs[ssum>0])^.3,col="blue")
text(xyr[rsum>0,1],xyr[rsum>0,2],actors[rsum>0],cex=2.7*(rr[rsum>0])^.3,col="tan")

#text(xys[ssum>0,1],xys[ssum>0,2],actors[ssum>0],cex=1,col="blue")
#text(xyr[rsum>0,1],xyr[rsum>0,2],actors[rsum>0],cex=1,col="red")
